# This file contains the S-plus code that was used to implement the estimators and functionals in "Transformations in hazard rate estimation", Journal of Nonparametric Statistics, (20): 721-738, (2008)
# As a fully working example the code herein can reproduce Figure 4 of the paper (the relevant data are also given below). 
# The same functions (with different input data) can be used in reproducing other figures of the same paper 
# More general use of the estimates is also feasible by supplying the corresponding data set.  Use of the functions is self - explanatory and is further facilitated by the comments at he preamble of each function.

#N.B. Compatibility with R is not guaranteed as the functions were not written for use with R - however, with minimal modifications, their use in the R environment should be feasible.

x1<- c(1, 1, 1, 5, 7, 8, 8, 13, 14, 14, 17, 18, 21, 21, 22, 25, 27, 27, 30, 30, 31, 31,
32, 34, 35, 36, 37, 38, 39, 39, 40, 49, 49, 54, 56, 56, 62, 63, 65, 65, 67, 75, 76, 79, 82, 
83, 84, 84, 84, 90, 91, 92, 93, 93,103, 103, 111, 112, 119, 122, 123, 126, 129, 134, 144, 
147, 153, 163, 167, 175, 228, 231, 235, 242, 256, 256, 257, 311, 314, 322, 369, 415,  573,
 609, 640, 737)


"Epanechnikov"<-
function(x){
	#Epanechikov kernel, square root of 5 is taken as 2.236068 (to save comp. time).
	ifelse(abs(x)< 2.236068, (3/4)*(1-((x^2)/5 ))/2.236068, ifelse(abs(x)>2.236068, 0,0))
}

"IntEpanechnikov"<-function(x)
{
	ifelse(x< -2.236068, 0, ifelse(x> 2.236068, 1, .5- (2.236068*x*(x^2-15))/100))
}

"HigherOrder"<-function(x){ifelse(x < -4, 0, ifelse(x> 4, 0, dnorm(x, 0, 1) * (3-x^2)/2))}#, ifelse(abs(x) > 4, 0, 0)) }


"Thre1"<- #Watson - Leadbetter hazard rate (see Phd)
function(xin, xout, kfun, bandfactor)
{
    n<- length(xin)
    n1<-1:n
    xin.use<-sort(xin)
    h<-bandfactor*diff(quantile(xin.use, c(.25, .75)))#NOTE: It doesn't matter if I use xin.use or xin in the interquartile.
    arg1<-(sapply(xout, "-", xin.use))/h	
    arg2<-kfun(arg1)* 1/(n-n1+1)
    arg3<-colSums(arg2)* 1/h
#cat("######################################", h)
    arg3
}

"iThre1"<- #Watson - Leadbetter integrated hazard rate (cum. haz. est.)
function(xin, xout, kfun, bandfactor)
{
    n<- length(xin)
    n1<-1:n
    xin.use<-sort(xin)
    h<-bandfactor*diff(quantile(xin.use, c(.25, .75)))
    arg1<-(sapply(xout, "-", xin.use))/h	
    arg2<-kfun(arg1)* 1/(n-n1+1)
    arg3<-colSums(arg2)
    arg3 
}

"SimpsonInt"<- function(xin, h)
{   #xin = data, h=distance of the points at which function is evaluated
	 n<-length(xin)
    nn<-1:n
    int0<-xin[1]
    int2n<-xin[n]
    inteven<-xin[seq(3,n-1,by=2)]
    intodd<-xin[seq(2,n-2,by=2)]
    sumodd<-sum(intodd)
    sumeven<-sum(inteven)
    res<-(h/3)*(int0+4*sumodd+2*sumeven+int2n)
    res
}

"ikde"<- function(xin, xout, kfun, bandfactor)
{
    n<- length(xin)
    h<-bandfactor*diff(quantile(xin, c(.25, .75)))
    arg1<-(sapply(xout, "-", xin))/h	
    arg2<-kfun(arg1)
    arg3<-colSums(arg2)
    arg3/n
}

"Thre2"<-
	#Transformed hazard rate estimator (2nd stage) using nonparametric transformation -log(1-\hat F(x))
function(xin, xout, kfun, ikfun, bandfactor)
{
	n<-length(xin)
	n1<-1:n
	nn<-length(xout)
	xin.use<-sort(xin)
	hBand<-diff(quantile(xin, c(.25, .75)))
	GofX<- -logb(1- ikde(xin.use, xout, ikfun, bandfactor), base=exp(1))#, base=exp(1)) 	
	Tsample<- -logb(1- ikde(xin.use, xin.use, ikfun, bandfactor), base=exp(1))	
	h2<- .0013* hBand #.0019  .00135
	cat("44444444444444444444444444", h2)
	GofXdashed<- Thre1(xin.use, xout,  kfun, bandfactor)  				
	arg1<-(sapply(GofX, "-", Tsample))/h2
	arg2<-kfun(arg1)* 1/(n-n1+1)
	arg3<-colSums(arg2)
	arg4<-arg3*GofXdashed/h2 					
	arg4
}

x<-seq(1, 700,length=100)

arg1<-Thre2(x1, x, Epanechnikov, IntEpanechnikov, .39)

plot(x, arg1, type="l",  xlab = "x", ylab="Hazard rate")  # TRANSFORMED ESTIMATE HAS BEEN CREATED HERE
key(text=list(c("TKHRE", "VKHRE", "HKHRE"), font =4), 
lines = list(lty=c(0,3,6)), corner=c(0,0),
cex= 1.15) #cex specifies font size


"VarHre"<-
#Adaptive variable bandwidth hazard rate estimator
function(xin, xout, h2, bandfactor, kfun=Epanechnikov)
	{
	n<- length(xin)
	nn<-length(xout)
	n1<-1:n
	xin.use<-sort(xin)
	kernel<-Thre1(xin.use, xin.use, kfun, bandfactor) #Pilot estimate (W-L)
	SqrtKernel<-sqrt(kernel)
	vhre<-sapply(1:nn, function(i, xin.use, xout, h2, kfun, n, n1, SqrtKernel)
							{ 
							 sum( (SqrtKernel *  kfun( ((xout[i]-xin.use)/h2) * SqrtKernel)	)/(n-n1+1))							
							}, xin.use, xout, h2, kfun, n, n1, SqrtKernel)		
	vhre/h2
	}
	

lines(x, VarHre(x1, x, 1.8, .5),lty=3)#, xlab = "x", ylab="Hazard rate")  
lines(x, Thre1(x1, x, HigherOrder, .39), lty=6)